Absence of a Spin Liquid Phase in the Hubbard Model on the Honeycomb Lattice 
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A spin liquid is a novel quantum state of matter where a finite charge gap exists although the 
conventional band theory would predict metallic behavior. Finding a stable spin liquid in two or 
higher spatial dimensions is one of the most challenging and debated issues in condensed matter 
physics. Very recently, it has been reported that a model of graphene, i.e., the Ifubbard model on 
the honeycomb lattice, can show a spin liquid ground state in a wide region of the phase diagram, 
between a semi-metal and an antiferromagnetic insulating phase. Here, by performing numerically 
exact quantum Monte Carlo simulations, we extend the previous study to much larger clusters 
(containing up to 2592 sites), and find no evidence of this spin liquid region. Our calculations 
together with previously established results strongly indicate that it is possible to establish a spin 
liquid ground state only in different models with real frustration. 



A spin liquid is a Mott insulator that is not adiabat- 
ically connected to any band insulator. Recently much 
attention has been focused on a possible spin liquid in 
two or three spatial dimensions [iHal- On one hand, it 
has been suggested experimentally that several organic 
materials represent good candidates for spin liquids 
[^. On the other hand, it has been shown theoretically 
that only few very particular and simplified models in- 
deed display a spin liquid ground state [1, • 

In order to understand this issue, let us consider a 
model Hamiltonian on a lattice describing the insulat- 
ing state of electrons at half-filling, i.e., one electron per 
lattice site. Since the charge gap is assumed, only spin 
degrees of freedom remain and can be described by the 
spin-1/2 Heisenberg model. Since any spin-1/2 model 
corresponds to a well defined interacting hard core boson 
model, the crucial question is to understand how - at 
zero temperature - these bosonic degrees of freedom can 
avoid Bose-Einstein condensation, a necessary condition 
for a stable spin liquid with no long-range order of any 
kind. 

In this report, we study the ground state of the half- 
filled Hubbard model on the honeycomb lattice (see 
Fig. [1} defined by the following Hamiltonian, 

H ^ (cl^Cj^ +c]^c,^'j +uY^ni^n,i, (1) 

where {cia) is a creation (annihilation) operator of 
spin up/down (cr =t,i) electrons on lattice site z, Uia- — 
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FIG. 1: The honeycomb lattice. Primitive lattice vectors 
Ti and T2 are denoted by red arrows. As an example, the hon- 
eycomb lattice with L = 3 is indicated by dashed blue lines. 
Solid and open circles are on A and B sublattices, respectively. 



cl^Cia, and t (U) denotes the nearest-neighbor hopping 
(on-site repulsion). This is known to be a model Hamil- 
tonian for graphene with U/t « 3 0|. More impor- 
tantly, it is not geometrically frustrated, namely, as seen 
in Fig. [l] the neighboring sites of any site on A sublat- 
tice belong to B sublattice (and vice versa). Indeed, it is 
well known that the ground state becomes an antiferro- 
magnetic insulator (AFMI), i.e., classically Neel ordered, 
from a semi- metal (SM) with increasing ll/t [l^ • 

Very recently, Meng et al. [2| has reexamined the 
ground state phase diagram of this model and found a 
possible spin liquid phase in the range 3.5 < U/t < 4.3 
between SM and AFMI. Their finding is rather surpris- 
ing because it is widely believed that a stable spin liquid 
occurs most likely in frustrated quantum systems where 
strong quantum fluctuations destroy the long-range mag- 
netic order even at zero temperature [l|. Their study was 
particularly successful because, with the auxiliary field 
technique [13| , there is no sign problem in the correspond- 
ing quantum Monte Carlo simulations, and an accurate 
finite size scaling was possible by using numerically ex- 
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act results for clusters containing up to 648 sites. So 
far, their results represent the most important numerical 
evidence for a possible spin liquid ground state in a "real- 
istic" electronic model in two dimensions (2D), because, 
to our knowledge, only a particularly simplified quantum 
dimer model on the triangular lattice [9] and the Kitaev 
model [Toj, built ad hoc to have an exact solution, allow 
a spin liquid ground state in 2D. Furthermore, their re- 
sults were considered a clear violation of the "Murphy's 
Law": in a too simple model, not vexed by the "fermion 
sign problem" , nothing interesting can occur Q . 

Here, by performing simulations for much larger clus- 
ters containing up to 2592 sites, we show that Bose- 
Einstein condensation concomitantly occurs once the in- 
sulating behavior sets in, confirming the more conven- 
tional Hartree-Fock transition from SM to AFMI [T^ . 
Although our results agree with the previous study for 
the same clusters up to 648 sites, we have reached the 
opposite conclusion. This reminds us similar claims on 
spin liquid behaviors in different systems in 2D [l^, [l^ , 
which have been corrected later on by much larger clus- 
ter simulations, showing instead antiferromagnetic long- 
range order 



Results 

We use finite size clusters of iV = 2L^ sites (thus con- 
taining L X L unit cells) with periodic boundary con- 
ditions (see Fig. [T]), which satisfy all symmetries of the 
infinite lattice JJij (also see Supplementary information). 
Here L is the linear dimension of the cluster and we take 
L up to 36. We use the well established auxiliary field 
Monte Carlo technique 



13[, which allows the statistical 



evaluation of the following quantity. 



0(r) = 



(2) 



where O is a physical operator, \ipii) (|V'l)) is the right 
(left) trial wave function (not orthogonal to the exact 
ground state), and r is the projection time. The exact 
ground state expectation value (O) of the operator O 
is then obtained by adopting the limit of r — >■ oo and 
At — >■ for 0(r), where At is the short time discretiza- 
tion of T. This approximation is necessary to introduce 
the auxiliary fields and implies a systematic error, 
negligible for small At (see Supplementary information). 

First, we study both the spin structure factor S'af = 
{[w^r i'^r.A — Sr,B)] ) and the spin-spin correlations 
Cs{R} — {Sr.A - Sr+R.A) at the maximum distance \R\ = 
Lmax of each cluster for U/t = 4, where the strongest 
evidence of a spin liquid behavior was found in Ref. Q- 
Here Sr.A {Sr,B) is the spin operator at unit cell r on 
A (B) sublattice. As shown in Fig. [2)d, our results show 
consistently a finite value of the antiferromagnetic order 
parameter = Saf/N = C(Linax) for N oo, in 



sharp contrast to the existence of a spin liquid, i.e., spin 
disordered, ground state reported in Ref. 0. 

By doing similar calculations for several U/t values 
(see Fig. [5]), we find in Fig. [3] that approximately 
scales linearly with respect to U/t, i.e., nis oc \U — Uc\, 
which is similar to the critical behavior predicted by the 
Hartree-Fock theory . Corrections to this linear criti- 
cal behavior are obviously expected, but they should not 
change much the critical value to the antiferromagnetic 
order, which is estimated here to be Uc/t = 3.76 ± 0.04, 
much smaller than the one (~ 4.3) reported in Ref. S 

Let us now evaluate the spin gap Ag. In order to avoid 
possible errors in extrapolating the imaginary time dis- 
placed spin-spin correlation functions, here we calculate 
directly the total energies in the singlet and the triplet 
sectors, with improved estimators, which dramatically re- 
duce their statistical errors [l^ (also see Supplementary 
information). We can see clearly in Fig. |4^ that the ex- 
trapolated spin gaps for different U/t values are always 
zero within statistical errors (e.g., the statistical error as 
small as 0.004t for U/t = 4). 

Next, we investigate whether the system is metallic or 
insulating, namely, whether there exists a zero or a fi- 
nite charge gap. For this purpose, it is enough to simply 
study the long distance behavior of charge-charge corre- 
lations, p{R) = {nr.A'iT'r+R.A) ^ 1- Hcrc n^.A is the den- 
sity operator at unit cell r on A sublattice (see Fig[T]). 
They should change from power law to exponential be- 
havior at a critical U where the charge gap opens up. 
This change of behavior is evident in Fig. 0)3 in a C/ 
range of 3.7-3.8t, which is consistent with the onset of 
the antiferromagnetic transition (Uc), remarkably within 
a quite small uncertainty of < O.lt. Our results therefore 
strongly support the more conventional scenario of a di- 
rect and continuous quantum phase transition between 
SM and AFMI [IJ]. 



Discussion 

Let us now discuss here why we have not found any 
evidence of a spin liquid phase. As shown in Ref. [20, 
by applying one of the theorems by Lieb [llj, it is eas- 
ily proved that the exact ground state of this model for 
U ^ satisfies the Marshall sign rule [l^ in the sector 
of no doubly occupied sites, accounting for low energy 
spin excitations. Indeed, the phases coincide with those 
of the simple antiferromagnetic Neel state ordered along 
the x-spin quantization axis, which can be written as 

n (I t)R - I Dr) n (I t)R + I Dr), where | t)^ and 
ReA ReB 

I Dr are spin configurations (along the z-spin quanti- 
zation axis) at site R. The expansion of this state in 
terms of | t)jt and | Dr yields the simple Marshall sign, 
namely, it is negative if the number of spin down in the 
A sublattice is odd. Thus, the phases of the ground state 
are trivial in the bosonic spin 1/2 sector. Therefore, 
Bose-Einstein condensation can hardly be avoided and 
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FIG. 2: Finite size scaling of spin-spin correlation functions for the Hubbard model on the honeycomb lattice 
at half-filling. Spin structure factor, Saf, and spin-spin correlations at the maximum distance, Cs(Lniax), are denoted 
by triangles and circles, respectively. Here, L is the linear size of clusters containing A'' — 2L^ sites. Antiferromagnetic order 
parameter squared, m^, is estimated by finite size extrapolating Saf and Cs(-Lmax) to L — !• cxd, namely, ml — liniL^cx, Saf/N = 
limL_>oo Cs(Z/max). Solid curves are quadratic fit of the largest cluster data. It is clearly seen that a consistent extrapolated 
value rrig is obtained for both quantities Saf and Cs (Lmax), indicated respectively by triangles and circles at 1/L = 0. Error 
bars of the extrapolated values are computed with the bootstrap technique. Insets show the expanded plots for large L. 



a magnetic long-range order occurs once the charge gap 
becomes finite. 

At this point, one could be tempted to assume the gen- 
eral validity of the above observation for generic 8=1/ 2 
systems and use this criteria based on the phases of the 
ground state wave function as a powerful guideline in 
the search of spin liquids for model systems as well as 
for real materials. Indeed, in all unfrustrated spin- 1/2 
Heisenberg and Hubbard models in the sector of no dou- 
bly occupied sites, the phases of the ground state wave 
function are not at all entangled in real space as they 
factorize into independent contributions relative to each 
site. Instead, the phases of the ground state wave func- 
tions are highly non trivial in well established spin liq- 
uid models such as, for instance, the Kitaev's model [10| . 



and the celebrated quantum dimer model on the triangu- 
lar lattice , because they are described by paired wave 
functions, which couple in a non trivial way the phases 
of nearest neighbor spins [23j . 



Therefore, we conclude that in a true spin liquid in 2D, 
the phases of the ground state wave function should be 
highly non trivial and entangled, otherwise Bose-Einstein 
condensation would very likely destabilize any seed of 
spin liquid behavior. To our best knowledge, the above 
observation is valid so far for all spin-1/2 models, as we 
have ruled out the present case of the Hubbard model on 
the honeycomb lattice, and we do not know any excep- 
tion, not even from the experimental point of view. 
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FIG. 3: The ground state phase diagram for the 
half-filled Hubbard model on the honeycomb lat- 
tice. Antiferromagnetic order parameter rUs as a function 
of U /t. rris (indicated by solid circles) is obtained by finite- 
size extrapolating the root mean square of Saf/N, rUs = 
limi^oo Saf/N, as shown in Fig. [2] SM and AFMI stand 
for semi-metal and antiferromagnetic insulator, respectively. 
Solid line is the linear fit of data for = 4, 4.3, and 4.6. 



As = E{S = 1) — E{S = 0) with high accuracy, without 
facing the negative sign problem, by directly simulating 
the singlet 5 = and the triplet S = 1 sectors separately 
(see Supplementary information). 

In order to accelerate the convergence to the ground 
state, we use for \^Pr) sl Slater determinant with a def- 
inite spin S, by breaking only spatial symmetries to re- 
move the degeneracy at momenta K and K' for the clus- 
ters chosen (a similar strategy was adopted in Ref. Q). 
Conversely, we use for IiJjl) a rotational and translation- 
ally invariant Slater determinant obtained by diagonal- 
izing a mean field Hamiltonian, containing an explicit 
antiferromagnetic order parameter directed along the x- 
quantization axis. In this way, the left and the right trial 
wave functions break different symmetries (spin and spa- 
tial ones, respectively), and for any symmetric operator 
O the convergence to the ground state is expected to be 
much faster because it is dominated by the much larger 
gap Agap in the symmetric sector. Since Agap is expected 
to scale to zero (if indeed zero) at most as ~ l/L, we use 
a projection time t = {L + 4)/f, which we have tested 
carefully to be large enough for well converged results 
(see Supplementary information). We have also checked 
that the systematic error due to discretizing r is basically 
negligible with Art — 0.14 for spin gap calculations and 
with Art = 0.1 for correlation functions. 



Methods 



Here we employ the standard auxiliary field Monte 
Carlo algorithm [13| with a more efficient implementa- 
tion [l^ by using different left and right trial functions 
Ii/jl) and \tpR) in equation This also allows to in- 
clude in the trial wave function a Gutzwiller type corre- 
lation, exp(— (7 ni^Tiii), where g is the Gutzwiller vari- 
ational parameter, to optimize the efficiency. As reported 
in Ref. |l9l the statistical error in evaluating the energy 
E{S) for a given spin S is dramatically reduced for ap- 
propriate values of g. Thus we can evaluate the spin gap 
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scales to zero with increasing the resolution in momentum space, namely as 1/L. In the antiferromagnetic region, the spin 
gap should instead vanish as 1/L^. This explains why for U/t = 4.3 the gap extrapolates to negative values, as we are well 
inside the antiferromagnetic phase (see Fig. In any case, a sizable spin gap is not found for any value of U/t. b. Charge- 
charge correlation function p{R) = {nr,Anr+R,A) — 1 at the maximum distance \R\ = Lmax for several values of U/t. In the 
semi-metallic phase, p{R) ~ and L*' p{Lmsx) should converge to a finite value for L — >■ oo. Instead, when a charge gap 

opens, the charge-charge correlations should decay exponentially and L'^ p{Lmsx) converges to zero in this limit. A change of 
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